#include "curve_complex.h"

double func_ellipse(double a, double La, double Lb)
{//ellipse
  return (La*Lb/sqrt(Lb*Lb*cos(a)*cos(a)+La*La*sin(a)*sin(a)));
}

double func_ellipse_f(double a, double La, double f)
{//ellipse with origin at the focus, which is at the left of the center.
 //La: semi axis along long side.
 //f: distance from the center to the focus
 return (La*La - f*f)/(La-f*cos(a)); //focus at the left of the center
 //return (La*La - f*f)/(La + f*cos(a)); //focus at the right of the center
  
}
double func_circle(double a, double r, double c1)
{//circle
  return r;
}
double func_parabolic_f(double a, double f, double c1)
{//parabola with focus point at origin. and opening toward left.
 //a: angle in radian
 //f: distance from vertex to the focus.
  return 2*f/(1+cos(a));
}
double func_hyperbolic_f(double a, double La, double f)
{//hyperbolic curve (left side) with focus at origin and opening toward left.
 //a: angle in radian
 //La: vertex to the center
 //f: focus to the center
 double ecc = f/La;
 return La*(ecc*ecc-1)/(1+ecc*cos(a));
}
double func_spiral_Archi(double a, double Ca, double Cb)
{//general Archimedean spiral: r=Ca + Cb * a
//a: angle in radian
//Ca, Cb: coefficient
  return (Ca + Cb*fabs(a));
}
double func_spiral_log(double a, double Ca, double Cb)
{//general logarithmic spiral r=Ca*exp(Cb*theta)
 //a: angle in radian
  return (Ca*exp(Cb*a));
}

//support to call a function with a form (double func(double) ) to define the curve.
void arc_func(FUNC func, double a0, double da, double w, double c0, double c1)
{//support any function (double func(double, double, double )) defined curve
  const double eps = getlayer_accuracy(currentlayer());
  double r;
  double A, dA, Astop;  
  int ii;
    
  if(da<0)
  {
    da = fabs(da);
    a0 -=da;
  }
  A = a0;
  Astop = A + da;
  push_cp();//cp at the center/originate.
  
  ii = 0;
  mskpolystart(clayer, w);

  while(A<Astop)
  {
    //radius/curvature
    r = func(A, c0, c1);
    org(r*cos(A), r*sin(A), 0.0, NOFLIP);
    mskpolypoint(clayer,cp.x,cp.y);
    ii += 1;
    
    if (ii > (MAXPTS-1))
    {
      mskpolyend(clayer);
      
      //new polygon
      mskpolystart(clayer, w);
      mskpolypoint(clayer,cp.x,cp.y);
      ii = 1;//initiate the counting
    }
    //angle step
    dA = 2 * acos(r/(r + eps));
    A += dA;
    //fprintf(stderr, "%d, %e, %e, %e, %e\n", ii, A, dA, r, eps);   
  }
  //insert the last point (La*cos(Astop), Lb*sin(Astop))
  A = Astop;
  r = func(A, c0, c1);
  org(r*cos(A), r*sin(A), 0.0, NOFLIP);
  mskpolypoint(clayer,cp.x,cp.y);
  
  mskpolyend(clayer);
  pop_cp();  //final cp stays at the origin.  
}

void pie_func(FUNC func, double a0, double da, double c0, double c1, double x0, double y0)
{//pie surrounded by the arc defined by the func and the point (x0, y0)
  const double eps = getlayer_accuracy(currentlayer());
  double r;
  double A, dA, Astop;
  double x1, y1, xi, yi;
  
  int ii;
  bool bSplit = false;
  
  
  if(da<0)
  {    
    da = fabs(da);
    a0 -=da;
  }
  da = da - int(da/(2*M_PI))*2*M_PI;
  if(da<VERYSMALL) da= 2*M_PI;
  
  A = a0;
  Astop = A + da;
  push_cp();//cp at the center/originate.
  
  r = func(A, c0, c1);
  x1 = r*cos(A);
  y1 = r*sin(A);
  ii = 0;
  mskpolygonstart(clayer);

  while(A<Astop)
  {
    //radius/curvature
    r = func(A, c0, c1);  
    xi = r*cos(A);
    yi = r*sin(A);
    org(xi, yi, 0.0, NOFLIP);
    mskpolypoint(clayer,cp.x,cp.y);
    ii += 1;
    
    if (ii > (MAXPTS-3))
    {
      //insert the origin (0, 0) to close the arc 
      org(x0, y0, 0.0, NOFLIP);
      mskpolypoint(clayer,cp.x,cp.y);
      //insert the starting point
      org(x1, y1, 0.0, NOFLIP);
      mskpolypoint(clayer,cp.x,cp.y);
      //end the polygon
      mskpolygonend(clayer);
      //set the flag for the splitting
      bSplit = true;
      
      //new polygon
      mskpolygonstart(clayer);
      x1 = xi;
      y1 = yi;
      // use the last point as the starting point for new polygon
      org(x1, y1, 0.0, NOFLIP);
      mskpolypoint(clayer,cp.x,cp.y);
      ii = 1;//initiate the counting
    }
    //angle step
    dA = 2 * acos(r/(r + eps));
    A += dA;
    //fprintf(stderr, "%d, %e, %e, %e, %e\n", ii, A, dA, r, eps);   
  }
  //insert the last point (La*cos(Astop), Lb*sin(Astop)) to close the ellipse
  A = Astop;
  r = func(A, c0, c1);  
  org(r*cos(A), r*sin(A), 0.0, NOFLIP);
  mskpolypoint(clayer,cp.x,cp.y);
  
  //insert the origin (0, 0) to close the arc 
  org(x0, y0, 0.0, NOFLIP);
  mskpolypoint(clayer,cp.x,cp.y);
  //insert the starting point
  org(x1, y1, 0.0, NOFLIP);
  mskpolypoint(clayer,cp.x,cp.y); 
    
  //end the polygon
  mskpolygonend(clayer);
  pop_cp();  //final cp stays at the origin.
}

void ellipse(double La, double Lb, double w)
{
  if(w<=0)
    //pie_ellipse(La, Lb, 0.0, 2*M_PI);
    pie_func(func_ellipse, 0, 2*M_PI, La, Lb, 0.0, 0);
  else
    //arc_ellipse(La, Lb, 0.0, 2.01*M_PI, w);//introduce ~2C overlap to avoid weidge slot.
    arc_func(func_ellipse, 0, 2.01*M_PI, w, La, Lb);
  
//La: semi axis along x; Lb: semi axis along y;
 /* const double eps = getlayer_accuracy(currentlayer());
  double r;
  double A, dA;
  double x1, y1, xi, yi;
  
  int ii;
  bool bSplit = false;
    
  A = 0;
  push_cp();//cp at the center/originate.
  
  x1 = La;
  y1 = 0;
  ii = 0;
  mskpolygonstart(clayer);

  while(A<2*M_PI)
  {
    xi = La*cos(A);
    yi = Lb*sin(A);
    org(xi, yi, 0.0, NOFLIP);
    mskpolypoint(clayer,cp.x,cp.y);
    ii += 1;
    
    if (ii > (MAXPTS-3))
    {
      //insert the origin (0, 0) to close the arc 
      org(0, 0, 0.0, NOFLIP);
      mskpolypoint(clayer,cp.x,cp.y);
      //insert the starting point
      org(x1, y1, 0.0, NOFLIP);
      mskpolypoint(clayer,cp.x,cp.y);
      //end the polygon
      mskpolygonend(clayer);
      //set the flag for the splitting
      bSplit = true;
      
      //new polygon
      mskpolygonstart(clayer);
      x1 = xi;
      y1 = yi;
      // use the last point as the starting point for new polygon
      org(x1, y1, 0.0, NOFLIP);
      mskpolypoint(clayer,cp.x,cp.y);
      ii = 1;//initiate the counting
    }
    //radius/curvature
    r = fabs(pow(Lb*Lb*cos(A)*cos(A)+La*La*sin(A)*sin(A), 1.5)/(La*Lb));
    //angle step
    dA = 2 * acos(r/(r + eps));
    A += dA;
    //fprintf(stderr, "%d, %e, %e, %e, %e\n", ii, A, dA, r, eps);   
  }
  //insert (La, 0) to close the ellipse
  org(La, 0, 0.0, NOFLIP);
  mskpolypoint(clayer,cp.x,cp.y);
  
  if(bSplit)
  {
    //insert the origin (0, 0) to close the arc 
    org(0, 0, 0.0, NOFLIP);
    mskpolypoint(clayer,cp.x,cp.y);
    //insert the starting point
    org(x1, y1, 0.0, NOFLIP);
    mskpolypoint(clayer,cp.x,cp.y); 
  }
  //end the polygon
  mskpolygonend(clayer);
  pop_cp();  //final cp stays at the origin.
  */
}

void pie_ellipse(double La, double Lb, double a0, double da)
{//draw a pie based on ellipse edge with origin at center (polar angle)
 //La: semi axis along x; Lb: semi axis along y;
 //a0: initial pole angle
 //da: rotate pole angle, positive->anti-clocwise
 
 pie_func(func_ellipse, a0, da, La, Lb);
}

void pie_ellipse_f(double La, double f, double a0, double da)
{//draw a pie based on ellipse edge with origin at foucs (polar angle).
 //La: semi axis along x; Lb: semi axis along y;
 //a0: initial pole angle
 //da: rotate pole angle, positive->anti-clocwise
 
 pie_func(func_ellipse_f, a0, da, La, f);
}

void pie_parabolic_f(double f, double a0, double da)
{
  pie_func(func_parabolic_f, a0, da, f);
}
void pie_hyperbolic_f(double La, double f, double a0, double da)
{
  pie_func(func_hyperbolic_f, a0, da, La, f);
}

void arc_ellipse(double La, double Lb, double a0, double da, double w)
{//draw an arc based on ellipse edge  with origin at center (polar angle) 
 //La: semi axis along x; Lb: semi axis along y;
 //a0: initial polar angle
 //a: rotate polar angle
 //w: line width
 arc_func(func_ellipse, a0, da, w, La, Lb);
}

void arc_ellipse_f(double La, double f, double a0, double da, double w)
{//draw an arc based on ellipse edge  with origin at a focus point (polar angle)
 //La: semi axis along x; Lb: semi axis along y;
 //a0: initial polar angle
 //a: rotate polar angle
 //w: line width
 arc_func(func_ellipse_f, a0, da, w, La, f);
}

void arc_parabolic_f(double f, double a0, double da, double w)
{
  arc_func(func_parabolic_f, a0, da, w, f);
}
void arc_hyperbolic_f(double La, double f, double a0, double da, double w)
{
  arc_func(func_hyperbolic_f, a0, da, w, La, f);
}
void arc_spiral_Archi(double Ca, double Cb, double a0, double da, double w)
{
  arc_func(func_spiral_Archi, a0, da, w, Ca, Cb);
}

/*commented on 20121104, replaced by a simpler version
void pie_ellipse(double La, double Lb, double a0, double a)
{//draw a pie based on ellipse edge with origin at center (polar angle)
 //La: semi axis along x; Lb: semi axis along y;
 //a0: initial pole angle
 //a: rotate pole angle
 //
 
  const double eps = getlayer_accuracy(currentlayer());
  double r;
  double A, dA, Astop;
  double x1, y1, xi, yi;
  double LaLb, La2, Lb2;
  
  int ii;
  bool bSplit = false;
  LaLb = fabs(La*Lb);
  La2 = La*La;
  Lb2 = Lb*Lb;
  
  
  if(a<0)
  {    
    a = fabs(a);
    a0 -=a;
  }
  a = a - int(a/(2*M_PI))*2*M_PI;
  if(a<VERYSMALL) a= 2*M_PI;
  
  A = a0;
  Astop = A + a;
  push_cp();//cp at the center/originate.
  
  r = LaLb/sqrt(Lb2*cos(A)*cos(A)+La2*sin(A)*sin(A));
  x1 = r*cos(A);
  y1 = r*sin(A);
  ii = 0;
  mskpolygonstart(clayer);

  while(A<Astop)
  {
    //radius/curvature
    r = LaLb/sqrt(Lb2*cos(A)*cos(A)+La2*sin(A)*sin(A));  
    xi = r*cos(A);
    yi = r*sin(A);
    org(xi, yi, 0.0, NOFLIP);
    mskpolypoint(clayer,cp.x,cp.y);
    ii += 1;
    
    if (ii > (MAXPTS-3))
    {
      //insert the origin (0, 0) to close the arc 
      org(0, 0, 0.0, NOFLIP);
      mskpolypoint(clayer,cp.x,cp.y);
      //insert the starting point
      org(x1, y1, 0.0, NOFLIP);
      mskpolypoint(clayer,cp.x,cp.y);
      //end the polygon
      mskpolygonend(clayer);
      //set the flag for the splitting
      bSplit = true;
      
      //new polygon
      mskpolygonstart(clayer);
      x1 = xi;
      y1 = yi;
      // use the last point as the starting point for new polygon
      org(x1, y1, 0.0, NOFLIP);
      mskpolypoint(clayer,cp.x,cp.y);
      ii = 1;//initiate the counting
    }
    //angle step
    dA = 2 * acos(r/(r + eps));
    A += dA;
    //fprintf(stderr, "%d, %e, %e, %e, %e\n", ii, A, dA, r, eps);   
  }
  //insert the last point (La*cos(Astop), Lb*sin(Astop)) to close the ellipse
  A = Astop;
  r = LaLb/sqrt(Lb2*cos(A)*cos(A)+La2*sin(A)*sin(A));  
  org(r*cos(A), r*sin(A), 0.0, NOFLIP);
  mskpolypoint(clayer,cp.x,cp.y);
  
  if(bSplit || fabs(a - 2*M_PI)>VERYSMALL)
  {
    //insert the origin (0, 0) to close the arc 
    org(0, 0, 0.0, NOFLIP);
    mskpolypoint(clayer,cp.x,cp.y);
    //insert the starting point
    org(x1, y1, 0.0, NOFLIP);
    mskpolypoint(clayer,cp.x,cp.y); 
  }
  //end the polygon
  mskpolygonend(clayer);
  pop_cp();  //final cp stays at the origin.
}
*/

/*commented on 20121104, replaced by a simpler version
void pie_ellipse2(double La, double Lb, double a0, double a)
{//draw a pie based on ellipse edge  with origin at center (parametrical angle)
 //La: semi axis along x; Lb: semi axis along y;
 //a0: initial parametric angle
 //a: rotate parametric angle
 //
  const double eps = getlayer_accuracy(currentlayer());
  double r;
  double A, dA, Astop;
  double x1, y1, xi, yi;
  
  int ii;
  bool bSplit = false;
  
  
  if(a<0)
  {
    a0 -=a;
    a = fabs(a);
  }
  a = a - int(a/(2*M_PI))*2*M_PI;
  if(a<VERYSMALL) a= 2*M_PI;
  
  A = a0;
  Astop = A + a;
  push_cp();//cp at the center/originate.
  
  x1 = La*cos(A);
  y1 = Lb*sin(A);
  ii = 0;
  mskpolygonstart(clayer);

  while(A<Astop)
  {
    xi = La*cos(A);
    yi = Lb*sin(A);
    org(xi, yi, 0.0, NOFLIP);
    mskpolypoint(clayer,cp.x,cp.y);
    ii += 1;
    
    if (ii > (MAXPTS-3))
    {
      //insert the origin (0, 0) to close the arc 
      org(0, 0, 0.0, NOFLIP);
      mskpolypoint(clayer,cp.x,cp.y);
      //insert the starting point
      org(x1, y1, 0.0, NOFLIP);
      mskpolypoint(clayer,cp.x,cp.y);
      //end the polygon
      mskpolygonend(clayer);
      //set the flag for the splitting
      bSplit = true;
      
      //new polygon
      mskpolygonstart(clayer);
      x1 = xi;
      y1 = yi;
      // use the last point as the starting point for new polygon
      org(x1, y1, 0.0, NOFLIP);
      mskpolypoint(clayer,cp.x,cp.y);
      ii = 1;//initiate the counting
    }
    //radius/curvature
    r = fabs(pow(Lb*Lb*cos(A)*cos(A)+La*La*sin(A)*sin(A), 1.5)/(La*Lb));
    //angle step
    dA = 2 * acos(r/(r + eps));
    A += dA;
    //fprintf(stderr, "%d, %e, %e, %e, %e\n", ii, A, dA, r, eps);   
  }
  //insert the last point (La*cos(Astop), Lb*sin(Astop)) to close the ellipse
  org(La*cos(Astop), Lb*sin(Astop), 0.0, NOFLIP);
  mskpolypoint(clayer,cp.x,cp.y);
  
  if(bSplit || fabs(a - 2*M_PI)>VERYSMALL)
  {
    //insert the origin (0, 0) to close the arc 
    org(0, 0, 0.0, NOFLIP);
    mskpolypoint(clayer,cp.x,cp.y);
    //insert the starting point
    org(x1, y1, 0.0, NOFLIP);
    mskpolypoint(clayer,cp.x,cp.y); 
  }
  //end the polygon
  mskpolygonend(clayer);
  pop_cp();  //final cp stays at the origin.
}
*/

/* commented on 121104, replaced by a simpler version
void arc_ellipse(double La, double Lb, double a0, double a, double w)
{//draw an arc based on ellipse edge  with origin at center (polar angle) 
 //La: semi axis along x; Lb: semi axis along y;
 //a0: initial polar angle
 //a: rotate polar angle
 //w: line width
  const double eps = getlayer_accuracy(currentlayer());
  double r;
  double A, dA, Astop;  
  int ii;
  
  double LaLb, La2, Lb2;
  LaLb = fabs(La*Lb);
  La2 = La*La;
  Lb2 = Lb*Lb;
    
  if(a<0)
  {
    a = fabs(a);
    a0 -=a;
  }
  A = a0;
  Astop = A + a;
  push_cp();//cp at the center/originate.
  
  ii = 0;
  mskpolystart(clayer, w);

  while(A<Astop)
  {
    //radius/curvature
    r = LaLb/sqrt(Lb2*cos(A)*cos(A)+La2*sin(A)*sin(A));
    org(r*cos(A), r*sin(A), 0.0, NOFLIP);
    mskpolypoint(clayer,cp.x,cp.y);
    ii += 1;
    
    if (ii > (MAXPTS-1))
    {
      mskpolyend(clayer);
      
      //new polygon
      mskpolystart(clayer, w);
      mskpolypoint(clayer,cp.x,cp.y);
      ii = 1;//initiate the counting
    }
    //angle step
    dA = 2 * acos(r/(r + eps));
    A += dA;
    //fprintf(stderr, "%d, %e, %e, %e, %e\n", ii, A, dA, r, eps);   
  }
  //insert the last point (La*cos(Astop), Lb*sin(Astop))
  A = Astop;
  r = LaLb/sqrt(Lb2*cos(A)*cos(A)+La2*sin(A)*sin(A));
  org(r*cos(A), r*sin(A), 0.0, NOFLIP);
  mskpolypoint(clayer,cp.x,cp.y);
  
  mskpolyend(clayer);
  pop_cp();  //final cp stays at the origin.
}

void arc_ellipse2(double La, double Lb, double a0, double a, double w)
{//draw an arc based on ellipse edge  with origin at center (parametrical angle)
 //La: semi axis along x; Lb: semi axis along y;
 //a0: initial parametric angle
 //a: rotate parametric angle
 //
  const double eps = getlayer_accuracy(currentlayer());
  double r;
  double A, dA, Astop;  
  int ii;
    
  if(a<0)
  {    
    a = fabs(a);
    a0 -=a;
  }
  A = a0;
  Astop = A + a;
  push_cp();//cp at the center/originate.
  
  ii = 0;
  mskpolystart(clayer, w);

  while(A<Astop)
  {
    org(La*cos(A), Lb*sin(A), 0.0, NOFLIP);
    mskpolypoint(clayer,cp.x,cp.y);
    ii += 1;
    
    if (ii > (MAXPTS-1))
    {
      mskpolyend(clayer);
      
      //new polygon
      mskpolystart(clayer, w);
      mskpolypoint(clayer,cp.x,cp.y);
      ii = 1;//initiate the counting
    }
    //radius/curvature
    r = fabs(pow(Lb*Lb*cos(A)*cos(A)+La*La*sin(A)*sin(A), 1.5)/(La*Lb));
    //angle step
    dA = 2 * acos(r/(r + eps));
    A += dA;
    //fprintf(stderr, "%d, %e, %e, %e, %e\n", ii, A, dA, r, eps);   
  }
  //insert the last point (La*cos(Astop), Lb*sin(Astop))
  org(La*cos(Astop), Lb*sin(Astop), 0.0, NOFLIP);
  mskpolypoint(clayer,cp.x,cp.y);
  
  mskpolyend(clayer);
  pop_cp();  //final cp stays at the origin.
}
*/


/*commented on 121104, replaced by a simpler version
void arc_ellipse_f(double La, double Lb, double a0, double a, double w, bool f2AtRight)
{//draw an arc based on ellipse edge  with origin at a focus point (polar angle)
 //La: semi axis along x; Lb: semi axis along y;
 //a0: initial polar angle
 //a: rotate polar angle
 //w: line width
 //f2AtRight: if true, the second focus point is at the right side; otherwise, at the left side. 
  const double eps = getlayer_accuracy(currentlayer());
  double r;
  double A, dA, Astop;  
  int ii;
  
  double f = sqrt(La*La - Lb*Lb);
  double ecc = f/La;  
  double nom = La*(1-ecc*ecc);
  double sign = (f2AtRight)?-1:1;
    
  if(a<0)
  {    
    a = fabs(a);
    a0 -=a;
  }
  A = a0;
  Astop = A + a;
  push_cp();//cp at the center/originate.
  
  ii = 0;
  mskpolystart(clayer, w);

  while(A<Astop)
  {
    //radius/curvature
    r = nom/(1+sign*ecc*cos(A));
    org(r*cos(A), r*sin(A), 0.0, NOFLIP);
    mskpolypoint(clayer,cp.x,cp.y);
    ii += 1;
    
    if (ii > (MAXPTS-1))
    {
      mskpolyend(clayer);
      
      //new polygon
      mskpolystart(clayer, w);
      mskpolypoint(clayer,cp.x,cp.y);
      ii = 1;//initiate the counting
    }
    //angle step
    dA = 2 * acos(r/(r + eps));
    A += dA;
    //fprintf(stderr, "%d, %e, %e, %e, %e\n", ii, A, dA, r, eps);   
  }
  //insert the last point (La*cos(Astop), Lb*sin(Astop))
  A = Astop;
  r = nom/(1+sign*ecc*cos(A));
  org(r*cos(A), r*sin(A), 0.0, NOFLIP);
  mskpolypoint(clayer,cp.x,cp.y);
  
  mskpolyend(clayer);
  pop_cp();  //final cp stays at the origin.
}
*/